Installation

Install packages if not installed yet

# Install package. Enter "a" in console when prompted
if (!require("BiocManager", character.only = TRUE)) install.packages("BiocManager")
if (!require("EBImage", character.only = TRUE)) BiocManager::install("EBImage")
if (!require("tidyverse", character.only = TRUE)) install.packages("tidyverse")
if (!require("broom", character.only = TRUE)) install.packages("broom")

Load packages

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.1.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(broom)
library(EBImage)
## 
## Attaching package: 'EBImage'
## The following object is masked from 'package:purrr':
## 
##     transpose

Note

Image processing

Read a test image. I have cropped the data such that

# Read an image
img <- readImage(here::here("data/raw/test2.jpeg"))
display(img, method = "raster", all = T)

Reduce the dimension of the image from 3 to 1. Using the most contrast dimension

# Grey scale
colorMode(img) = Grayscale
display(img, method = "raster", all = T)

img <- img[,,3]
display(img, method = "raster")

Blur the image as the pre-treatment for gating out the background noise

# Blur
img_bl <- gblur(img, sigma = 1)
display(img_bl, method = "raster")

Set the global threshold

# Global threshold
threshold <- otsu(img_bl)
img_th <- combine(mapply(function(frame, th) frame > th, getFrames(img), threshold, SIMPLIFY=FALSE))
display(img_th, method = "raster")

Label the objects. Each object is a group of adjacent pixels.

img_label = bwlabel(img_th)

Gating out the background. Be cautious that here I am setting a manual threshold (e.g., objects number of pixel <= 10 are considered noise).

# Remove background
temp <- img_label %>% table()
index_rm <- names(temp)[which(temp <= 10)] # Pixel size for gating out background
img_clean <- rmObjects(img_label, index_rm)
display(img_clean, method = "raster")

Watershed separates the adjacent colonies

# Watershed
img_ws <- watershed(distmap(img_clean), 1)
display(colorLabels(img_ws), method = "raster")

You can get shape data for each colony: area, perimeter, radius, etc. The number of objects are the number of colonies.

colony_shape <- computeFeatures.shape(img_ws) %>% as_tibble()

Number of total colonies

nrow(colony_shape)
## [1] 34
colony_shape
## # A tibble: 34 × 6
##    s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
##     <dbl>       <dbl>         <dbl>       <dbl>        <dbl>        <dbl>
##  1    533         106         12.8         2.00         8.84         17.3
##  2    420         103         11.3         1.74         8.13         15.5
##  3    290          74          9.15        1.20         6.30         11.8
##  4    272          90          8.86        1.83         4.54         12.2
##  5    258          76          8.77        1.30         6.29         11.4
##  6    248          66          8.39        1.04         6.41         10.8
##  7    243          68          8.40        1.06         6.09         10.6
##  8    270          81          8.86        1.39         5.88         11.9
##  9    332          94         10.2         2.60         4.60         14.9
## 10    228          72          8.09        1.25         4.44         10.5
## # … with 24 more rows

Test on another image

# Read an image
img <- readImage(here::here("data/raw/test3.jpeg"))
display(img, method = "raster", all = T)

# Grey scale
colorMode(img) = Grayscale
display(img, method = "raster", all = T)

img <- img[,,3]
display(img, method = "raster")

# Blur
img_bl <- gblur(img, sigma = 1)
display(img_bl, method = "raster")

# Global threshold
threshold <- otsu(img_bl)
img_th <- combine(mapply(function(frame, th) frame > th, getFrames(img), threshold, SIMPLIFY=FALSE))
display(img_th, method = "raster")

img_label = bwlabel(img_th)
display(img_label, method = "raster")

# Remove background
temp <- img_label %>% table()
index_rm <- names(temp)[which(temp <= 10)] # Pixel size for gating out background
img_clean <- rmObjects(img_label, index_rm)
display(img_clean, method = "raster")

# Watershed
img_ws <- watershed(distmap(img_clean), 1)
display(colorLabels(img_ws), method = "raster")

In this case, the difference in colony size is very large. It makes sense to tell them apart by the size

colony_shape <- computeFeatures.shape(img_ws) %>% as_tibble()
colony_shape
## # A tibble: 50 × 6
##    s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
##     <dbl>       <dbl>         <dbl>       <dbl>        <dbl>        <dbl>
##  1   6069         368         44.1        3.73         37.4         53.5 
##  2   5467         317         41.4        2.79         35.4         49.4 
##  3   3517         242         33.6        4.76         24.8         41.9 
##  4   2616         218         28.8        4.07         18.5         38.7 
##  5   2664         254         28.9        1.31         26.0         32.5 
##  6    948         190         17.1        5.63          3.47        27.3 
##  7    261          66          8.70       0.961         7.10        11.4 
##  8    237          59          8.27       0.909         6.26        10.8 
##  9    236          60          8.27       0.813         6.19        10.2 
## 10    241          58          8.39       0.746         6.73         9.73
## # … with 40 more rows

Relative count by clustering

colony_shape %>%
    as.matrix() %>%
    kmeans(2) %>% 
    broom::tidy() %>%
    pull(size)
## [1] 45  5